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Abstract. We investigate the nature of the dynamically inactive phase of a 
simple symmetric exclusion process on a ring. We find that as the system's 
activity is tuned to a lower-than-average value the particles progressively lump 
into a single cluster, thereby forming a kink in the density profile. All dynamical 
regimes, and their finite size range of validity, are explicitly determined. 
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1. Introduction 

The study of temporal large deviations in simple systems such as those of 
interacting diffusing particles has revealed an unexpected wealth of behaviors, 
ranging from puzzling correspondences between driven systems and their equilibrium 
counterparts [1, 2, 3, 4], to intrinsically dynamic phase transitions [5, 6, 7, 8, 9, 10]. 
The dynamical phase transition terminology dates back to the eighties [11], in an 
era when the dynamical systems community exploited the so-called thermodynamic 
formalism of Ruelle and others to study and characterize the variety of dynamical 
regimes displayed by simple iterated maps. There have been many efforts in the 
recent past to classify the various phase transitions that can be found between distinct 
dynamical regimes, but there is a great scarcity of studies attempting an even coarse 
description of the properties of these dynamical regimes. Exceptions are Bodineau and 
Derrida's study of a weakly asymmetric exclusion process [5] , Toninelli and Bodineau's 
work on kinetically constrained models [12] and the finite size study of Bodineau, 
Lecomte and Toninelli [13]. In the generic dynamical phase transition scenario, one of 
the dynamical phases is disordered and its physics is easy to grasp, while the ordered 
symmetry-breaking phase is usually hard to characterize, for in dynamics there is no 
such thing as an easy-to-implement free-energy based variational principle. 

In the present work, we would like to devote our own efforts to the understanding 
of the nature of the various dynamical phases of the Simple Symmetric Exclusion 
Process (SSEP) on a ring, a system of mutually excluding particles hopping (if the 
target site is empty) with equal rates to either of their nearest neighbor sites on a 
one-dimensional lattice with periodic boundary conditions. The SSEP belongs to the 
broader class of diffusive systems whose dynamics can be described by fluctuating 
hydrodynamics, meaning that at the scale fixed by the system size, the field of 
occupation numbers becomes a smoothly varying function evolving through a Langevin 
equation with a vanishingly small noise (as the system size increases). We choose 
to classify the various time realizations the system can follow according to their 
activity level defined as the number of particle moves having taken place over a 
given time interval. In practice, however, we introduce a Lagrange multiplier s 
biasing trajectories towards a given average activity [14]. The s > (resp. s < 0, 
s ~ 0) regime allows one to explore the less-than-average-activity (resp. larger-than- 
average, typical) trajectories. The latter activity is a physical observable that has 
received considerable attention in the past [7, 14, 15], not only in diffusive systems 
but also in a variety of systems with slow dynamics [16, 17, 18], and which has strong 
connections with entropy production in noncquilibrium systems [19]. This is the 
simplest trajectory-dependent time-reversal symmetric observable one can think of. 
For realizations showing typical or larger than typical activity, the particles remain 
uniformly distributed, and the density profile remains flat. This is the disordered 
regime. By contrast, if we focus on realizations displaying a smaller than typical 
activity (s < 0), the uniform profile becomes unstable. This scenario was identified in 
[20] , and presents features similar to the transition at high enough values of the current 
in the weakly asymmetric exclusion process [6]. The transition point, occurring right 
at the typical (or average) activity is characterized by a universal scaling function that 
does not depend on the specifics of the SSEP, in the large size limit. The purpose 
of our work is precisely to investigate the properties of the ordered phase in which 
translation invariancc, as our intuition tells us, will be broken. The extreme low- 
activity behavior is indeed going to be dominated by realizations in which particles 
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have clustered into a single lump, thus permitting particle moves at the borders of the 
cluster, leading to a strongly subextensive activity. 

Before we proceed, we would like to recall that studying the statistics of the 
activity in the SSEP is equivalent to investigating the ground-state properties of a 
ferromagnetic XX Z chain [20], which allows one to connect our work to the integrable 
systems literature. We are indeed interested in the smallest eigenvalue —ip{s) of the 
operator 

H(.) - ^ - ^J21<tM+i + ^M+i + cV^al+J (1) 

fc=i 

where the erf's are the Pauli matrices. In standard notations (see Baxter [21] for 
instance) the anisotropy parameter A is related to the Lagrange multiplier s enforcing 
a given average activity by A = e''. Another related problem [22, 23] is that of 
determining the ground state of a gas of hard-core bosons with interaction strength 
— s and density po, which, in the present s > case, would thus correspond to 
attractive interactions (this is the no-saturation case discarded by Lieb and Liniger 
at the beginning of their seminal work [22]). In the regime of interest here, we have 
that A > 1. The same problem of attractive bosons with hardcore interactions is of 
direct interest to study the directed polymer in the replica approach [24], for which 
the problem was fully diagonalized [25, 26], see [27] for a review. It was shown [28] 
that such a model and such a regime for A were relevant for the description of 
the superconducting phase of some generalized Hubbard models, when working at 
fixed magnetization. The latter constraint, which renders the problem nontrivial 
(without the constraint, see section 8.8 in Baxter [21]) turns out to be our working 
framework, since particles are conserved in our closed ring, and magnetization per 
site is 2po — 1- While the active dynamical phase lent itself moderately easily to a 
Bethe Ansatz approach [20] , which was shown to be fully equivalent to the fluctuating 
hydrodynamics approach, it can be understood from [28] that things only get worse for 
A > 1. It therefore appears to be desirable to resort an alternative method. It seems to 
us that the phase separation mechanism conjectured in [28] has a simple interpretation 
in our approach, as will become clear later in the paper. Fluctuating hydrodynamics, 
in the integrable systems language, can be viewed as the correct effective field theory 
able to describe low lying excitations of the system. 

The questions we wish to answer can be phrased as follows. As the system's 
activity is lowered below its average, how does the flat density profile become unstable? 
Is there any limiting profile deep into this inactive phase? Which are the activity 
scales, and what are their system size dependence, governing these various dynamical 
inactive regimes? Given that these are generic questions that can be raised for other 
comparable systems, it is of interest to present a model with a dynamic phase transition 
for which they can be answered in an exact manner. We will begin with a brief 
overview of existing results related to large deviations in the SSEP, as predicted by 
fluctuating hydrodynamics — also termed Macroscopic Fluctuation Theory by Bertini 
et al. [29, 30, 31, 9, 10]. This will allow us, in Section 3, to present exphcit results on 
the density profile in the inactive regime. Three scaling regimes exist, which we will 
investigate separately in a mathematically controlled fashion. Open questions will be 
gathered in the final section. 
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2. The symmetric exclusion process and its fluctuating hydrodynamics 
description 

2.1. Fluctuating hydrodynamics 

Wc consider a Simple Symnietric Exclusion Process (SSEP) on an L-site one- 
dimensional lattice with periodic boundary conditions. Each of the mutually excluding 
N particles can hop with rate 1 to either of its nearest neighboring site, if empty. We 
consider time realizations of the process over the time interval [0, t] and we choose to 
work in space and time units scaled by the system size L and the diffusion relaxation 
time L^, respectively. We thus introduce a smoothly varying field p{x, r) of occupation 
numbers, defined by nj{t') = p{j/L,t' /L'^), whose existence we assume. The field p 
evolves according to a conserving Langevin equation 

drP + d,3 = (2) 

with j = —DdxP y=\/a£_, and where D{p) = 1 and a{p) = 2p{l — p) are functions 

reflecting the SSEP dynamics at a macroscopic scale. The Gaussian white noise, 
^(x, r), has unit correlations, (^(x, t)£^{x' ,t')) = 5{x — x')S{t — r'). The Green-Kubo 
relation ensures that — = /?/", where f{p) is the free energy per unit length, which, in 
the SSEP, reduces to the purely entropic expression — /?/ = —p\np—{l—p) ln(l— p) [up 
to p^, f" is the inverse isothermal compressibility, kt = (p^/")~^]- The conditions 
under which fluctuating hydrodynamics is valid, and the reasons why it applies to 
the SSEP, have been discussed by many authors. We refer the reader to Kipnis and 
Landim's book [32] which specifically addresses the SSEP, and to the series of papers by 
Bertini et al. [29, 30, 31, 9, 10], the latter employing the macroscopic fluctuation theory 
vocabulary. See also [2] for a field theory approach using coherent state construction 
of path integrals to represent the {c^'^^) = (e^*^(-'*)) [the operator H(s) is defined 
in (1)]. 

We will focus on the so-called activity K{t), which, in the SSEP, is the number of 
particle hops that have taken place over a given time window [0, t], over the whole ring. 
Up to irrelevant finite size corrections, we may view the activity as a functional of the 
local occupation numbers given by K{t) = J^ dt' J2j ± '^i(^')(l~'^i±i(*'))- ^^ terms of 

the field of occupation numbers, it is expressed by K{t) = L^ J^ dx J^ dTa{p{x, r)). 
The equilibrium distribution in the SSEP is a simple Bernoulli distribution with 
parameter po- Denoting by po = N/L the average density, we thus find that, to leading 
order in the system size, {K)/{tL) — 2 — \^i~^_i\ — <^{po)- For time realizations 
with activity greater than (K) we expect that the system remains as homogeneous 
as possible to favor particle hops able to accommodate a high value of K. We also 
expect, as we confirm below, that for low activities the system will group particles 
into a single cluster. We now briefly recall what the existing literature has established 
about these opposite regimes. 

2.2. Universal fluctuations of the activity 

The key to the calculation of large deviation properties in systems described by a 
Langevin equation (2) is that, since the noise becomes vanishingly small as the system 
size increases, a WKB-like saddle point expansion can be exploited. There arc many 
ways to implement the saddle point expansion. One is to cast the generating function 
Z{s,t,L,po) = (c^"^) of the activity in the form of a path-integral over a pair of 
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fields p and p (the latter is the actual occupation number field, while the former is 
conjugate to the noise £,). The net result of that procedure that has been described 
many times in the past (see 6.17. [2]) is that 

Sp&pc-^'^^P'P''^ (3) 



where the action S[p, p; s] reads 

pdrp + dxpdxP Tridxpf + sL?a{p) 



ft/L' rl r 

S^ At dx pdrP + d^pd^p-^{d^pf + sL'<j{p) (4) 

Jo Jo I ^ 

Once cast in the above path-integral formulation, it is clear that, as i — > 00, the 

leading behavior to Z will be given by a saddle point evaluation of the path-integral. 

One must then minimize the action S with respect to p and p. The latter gives us the 

optimal density profile able to accommodate a given value of the average activity, as 

tuned by s. In Appert et al. [20] it was shown that the large deviation function of the 

activity, defined by 

V'(s) = lim (5) 

t^toa t 

could be written in the form 

where J^ is a scaling function independent of the functions D and cr, and its argument 
reduces to J^(§sL^) in the SSEP. The functions D and a appearing in this formula 
are evaluated at poj a-s in the sequel of the paper when no argument is made explicit. 
The hypotheses underlying this result are that the optimal profile in the path-integral 
formulation for the generating function Z is both stationary and uniform at the value 
Po- The function T(u) exhibits a branch cut along the real axis for u -^ (7r^/2)~, 
which signals that for s > ^, Ac = ir'^ /(t{pq)^ the uniform profile ceases to be the one 
minimizing the action. The new results of the present work are devoted to a study of 
this sL^ > Ac regime. We now define the parameter A by A = sL^. 

It is piquant to note that T has different integral representations according 
to whether A < or < A < Ac. We need it on the positive side, namely 
J"(a;) = 2.T — 4-\/2a;^/^ /_^ dy j/^ cot (\/2x^/T^^y^) . The function J" has the following 
limiting behaviors 

97/2 
T{x) ~ -7—i-xf^^ for X -^ -00 (7) 

_2- 

T{x) = Ci - C2(7rV2 - a;)^/2 + , _ foj. x ^ — (8) 

where Ci ~ 51.61351... and C2 — 55.83091... > 0. The latter expansion explicitly 
displays the ^i{^ j^ — x singularity. The x — > — 00 limiting behavior was first found 
by Lieb and Liniger in their study of the one-dimensional Bose gas with repulsive 
interactions [22] . We now turn to our analysis of the inactive phase. 
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3. Inactive phase 

3.1. The technical problem 

Inspired by the recent work of Bodineau and Toninelh [12], we define a rescaled large 

deviation function ^(A) = linii_j.oo -^ linit_i.oo — — 2 = limi^.oo i^(A/i ). From 

(6) we get that 0(A) = — Acr(po) which has finite size corrections of the order 1/L given 
by F{a\/2). While F appears only through finite-size corrections, we note that its 
singularities signal the end of the regime over which the uniform saddle profile remains 
the optimal one. Mathematically, the quadratic expansion of the action around the 
uniform saddle ceases to be positive definite at A > Ac- We therefore search for another 
saddle to the Euler-Lagrange equations. Due to the conservation of the number of 
particles, we introduce a Lagrange multiplier ^ enforcing the total density to be poj 
by adding a fi J {p — po) contribution to the action. Assuming a stationary saddle, the 
Euler-Lagrange equations of motion read 

j = -dp + adp, 0*2/0 + Y {dp)^ - Ao-' - /i = (9) 

We see a posteriori that in the homogeneous regime we must have fi = — A(t'(po)- The 
only way to find out the stability regime of the homogeneous profile is to study the 
quadratic action expanded in the vicinity of that saddle. Since biasing the trajectories 
with the activity does not break the left-right symmetry nor the time-reversal one {K 
is even by time-reversal), we search for a solution to (9) in which the current j = 0. 
Then it is easy to verify that the equations in (9) are equivalent to Euler-Lagrange 

equations deduced from a Lagrangian C = ^g — h Act -I- (p — Po)Mi and were the 
momentum conjugate to p is tt = -mp- = ^^ — dxp- The corresponding Hamiltonian 
% = TrdxP — £ is a constant of motion, which means, when written in terms of the 
Lagrangian variables, that Ti, given by 

n^^{dpf-\a-{p-po)p (10) 

is independent of x. A simple rewriting of the above equation tells us that 

\{dxpf + Ep{p) = 0, Ep{p) = -\a{pf - [{p - po)fi + H]a{p) (11) 

Again, (11) has a simple mechanical interpretation: p is the position of a particle 
evolving in the force field given by the potential energy Ep{p), with kinetic energy 
i{dxp)^, while maintaining a zero total mechanical energy. Time coordinate is x. 
The two constants "H and p arc for the moment undetermined. Periodic boundary 
conditions and particle conservation will lead to A and po dependent expressions for 
these parameters. In the homogeneous regime (A < Ac), we must have Ep = 0, 
TL ~ ~Xa(po) and p ~ —Xa'{po). In the mechanical analogy, searching for a periodic 
profile means that we are after a periodically oscillating solution between two extreme 
positions. For A < Ac and the corresponding values of TL and p there is only a single 
value of p such that Ep{p) ~ 0, which means that the oscillations are reduced to 
a standstill at this very density (which is of course po). Note that, mathematically 
speaking, a similar set of equations were obtained by Bodineau and Derrida [33] in the 
context of the additivity principle for the large deviations of the current in a system 
in contact with two reservoirs (a situation in which particle number is not conserved 
and no chemical potential is required). We now set out to find the optimal profile p{x) 
along with H and p in the regime where A — >■ A+. 
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Figure 1. The potential energy profile, Ep{p), is shown for a set of typical values 
of the parameters A, H, fj, and po- 



3.2. Close to the critical point 

We use this formulation to find the optimal profile for A = Ac + e, with e ^ O"*". 
In the spirit of [5], we assume that a small perturbation around the flat profile 
will develop, which leads us to write that p{^x) = po + \/e4>i{^) + £02(a;) + . . ., 
IJ, = He + y/sHi + e/X2 + . . ., with fie = — tt^ — . We write that H is independent 
of X order by order in e, which, after a few manipulations, leads to 

Mx) = -[Po(l - Po)]^/' eos(2^2;) 



and 



MX 



-cr'(po)CT(po)^cos(47rx) 



together with /ii = and /i2 = —a-'{po)^ /A. We then find that as A = Ac 



0(A) = -7r2-a(po)e + 3 



o-(Po)' 



(12) 

(13) 
e, with 

(14) 



We find that H = — tt^ — cr(l — 4(T)e + . . .. In other words, the instability sets in via 
the largest wave-length Fourier mode, just as was the ease in the system studied by 
Bodineau and Derrida [5]. However, in our case, we can go deeper into the broken 
symmetry phase, as will now become clear. 



3.3. In the inactive regime 

In order to pursue our mechanical analogy, we plot Ep (p) as a function of p and have 
the space variable x play the role of a time. A typical plot is shown in figure 1. There 
one can see that Ep{p) vanishes at two densities p± given by 



P± = 



2X + p± yj%\[% - Hpo) + (2A + Hf 



4A 



(15) 
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The particle thus osciUates between these two extreme positions, and it takes a fuU 
period 



1 



1 fP+ 
dx = 1 = 2 / dp- , 

Jp_ V-'^Ep{p) 



(16) 



to do the round-trip between p_ and p_|-. Besides, we must have that the total density 
is po, which in turn imposes 



Po 



Axn(x) = 2 / dp — = 

7p_ ^V-2i?p(p) 



(17) 



The two equations (16) together with (17) thus completely fix the unknowns 'K and p. 
At this stage, it is more convenient to deal with the po = 1/2 case, for which analytic 
expressions are somewhat simpler. 



3.3.1. At half-filling Particle-hole symmetry at po = 1/2 imposes that Ep{p) = 
Ep(l — p), which in turn forces p = (this can explicitly be verified from (16) and 
(17)). Then we have that 



A - J\{2'H + A) A + J\{2n + A) ^ 
< /O- = ^^77 < P+ = ^'—T-. < 1 



(18) 



2A '^ 2A 

The density profile thus evolves between p_ and p+ given above, with the constant H 
being given by the implicit solution to 



2 " V^ 

where K is the complete elliptic integral of the first kind. Using that as w — )• 

1 u 

K(l-u) ==21n2- -lnu+-(-lnii-2 + 41n2)+0(u2) 

2 8 

we can solve for H as a function of A from (19), which leads to 
n = 8Ae"\/? -I- 16e"2^/?A(\/2A - 4) 

-t- 16e"3v^A(6A - 19\/2A + 22) + (x^^^c^^y^ 
The large deviation function 4'{X) is then obtained from 

"(a.p)2 



0(A) 



dx 



= / dx 
Jo 



n- 



2a- 
{d.P? 



+ Act + (p - po)p 



Given that we have 



o- Jp_ 



a 



we arrive, after some manipulations, at 



(/)(A) = -V2(V2-H + A + VA) E 



2y/X{mTX) 



(19) 
(20) 

(21) 

(22) 
(23) 

(24) 
(25) 



n + x + y'\i2n + \)J 

where "H as a function of A is extracted from the implicit formula (19). The notations 
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Xc + Yo ('''^'^i large dashes), 25 



Figure 2. The density profile is shown for A 
(blue), 200 (purple, small dashes) at average density po = o- '^^'^ kink becomes 
steeper as A is increased. The slope of the rise to the plateau increases as vA for 
large A. 



K and E refer to the complete elliptic integrals of the first and second kind, respectively. 
This expression is valid for any A > Ac- Using the following asymptotic formula for E 
as u — > 



E(l - u) = 1 



1, 1 

-Inu - - 
4 4 



In 2 



u 
64 



(-6 In It -13 + 24 In 2 
we arrive at the following large A asymptotics for 0(A): 
(/)(A) = - 2%/2A - 8\/2Ac"^/? 

+ 8(2A - 3\/2A)c"2Vl + o (xh 



0(m3) 



-3^?' 



(26) 



(27) 



This is the final result for the po = 1/2 particular case. Translated into the language of 
quantum spin chains, 0(A) in (25) is the ground-state energy per site of a ferromagnetic 
XX Z chain with anisotropy parameter A = 1 + -p- for A > Ac, and at half-filling. To 
our knowledge, this is the first appearance of this expression in the literature. The 
optimal profile corresponding to that value of 0(A) takes the form of a smooth kink 
with area equal to 1/2, as shown in figure 2 for increasing values of A. Its expression 
reads (for x as function of p for < a; < i) 



^/2F 



^ip) = 2 - 



(p-l)A-g 



il-p}y/X{2H+\} 



2y/X(2H+X} 



H+\+^\(2H+\) 



s/2H + A + VA 
where F is the incomplete elliptic integral of the first kind. 



(28) 
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3.3.2. At arbitrary density po For general poi by exactly the same methods, though 
with somewhat less elegant simplifications (no symmetry enforces the Lagrange 
multiplier /i to vanish), the two equations (16) and (17) now take the following forms 



1 
2^0 = 



2AK 



Ai+2A(l-2p) 
(,i+2A(l-p))(l-p) 



(,.-2Ap)n(H±|M 



Ai+2A(l~p) 



{p+2A(l-p))(l-p) 



where 11 is the complete elliptic integral of the third kind, along with 
' A* + 2A(1 - 2p) 



1 



K 



(1 - p)(2A(l -p)+^i) \ifi + 2A(1 - p))(l - p) 



(29) 



(30) 




Figure 3. Large deviation function (j>{X) (thick line) as a function of A, at average 



density po 



The function is singular at A = Ac 



O-(PO) 



(the second derivative 



0(A) = -2v/(l-/5)(2A(l-/5)+/i)E 



is discontinuous), marking the entrance in the regime where the density profile is 
non-uniform (A > Ac). The value of the large deviation function for a uniform 
profile at A > Ac is shown for comparison (dashed line). 

where p is either equal to the minimum or the maximum density p± [defined in (15)]. 
Once p and Ti have been extracted from (29) and (30), one can determine the large 
deviation function 0(A) as a function of A and po from the formula (see figures 3 and 4) 

p + 2X{l-2p) \ 
,(M + 2A(l-p))(l-p)y' ^ ' 
where A > Ac. The expansion close to Ac can be checked from those results. We see in 
particular that xk(i^^)c = 0"(A) presents a jump 0"(A+) - 0"(Aj) = ^^^^. We do 
not intend to labour beyond the leading order of (/)(A) as A ^ 1, which is independent 
of po, 

0(A)~-2V2A + -H + (1-po)m+--- (32) 

That the leading contribution is independent of po could have been expected, as the 
main contribution to (f> as given in (23) is the integral (24), which is itself dominated by 
a small region of size 1/vA where the derivative of p varies the most steeply between 
p_ and p+, while density is close to ^^ „^^ ~ 1/2, and thus 



0(A) 



i/\/a' 



l/VX 



Va 



(33) 
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0.4 



02 



0.1 




20 



Figure 4. Mean activity 7t(^)a = ~<t>' {^) (thiclc blue line) as a function of A, at 



average density po 



The function is singular at A = Ac 



-2-^ . The value of 

-lj{K)s is shown for a uniform profile at A > Ac is shown for comparison (dashed 
line). As expected, above Ac, the activity for histories with a non-uniform profile 
is lower than for uniform profiles. 



Together with the impUcit formulas (29) and (30), equatfon (31) is the ground 
state energy of an XXZ ferromagnetic Heisenberg chain with anisotropy parameter 
A = 1 + -p- at fixed magnetization L{2po — 1). The density profile is, for x as function 
of /9 for < cc < i (see figure 5), 



x{p) = 



arcsm 



in f I . (M-2Ap-)(p-p-) 



2A(l;-2p_2+£i 
(/^-2Ap_)p_ 



v/p-(2Ap_-M) 



(34) 



3-4-. The infinitely inactive regime 

In the infinitely inactive regime, that is for A 3> i^, the density profile described 
in the previous subsections becomes infinitely steep and the assumptions underlying 
fluctuating hydrodynamics fail. Since the typical slope of the rise of the profile towards 
its plateau value increases at vA, fluctuating hydrodynamics must fail when the slope 
exceeds L, that is when profile variations become significant at the lattice scale, which 
is consistent with the \':$> Li^ condition. For large s it is easy to see that the ground 
state of the s-modified evolution operator (1) reads 

4,{s) = -2La(po)(l - e-^) + 0{e-^') (35) 

The above regime consistently connects to the calculation of the previous subsection 
where it was shown that 0(A) ~ — a/A for large values of A. The validity of that 
previous result holds until 4>{X) ^ —2La{pa), when it connects with the infinitely 
inactive regime. This occurs as A ~^ i" and confirms once more that the crossover is 
at a typical A of the order of L^. 
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Figure 5. The density profile is siiown for A = Ac + ^ (red, large dashes), 25 
(blue), 200 (purple, small dashes) at average density Po = 3 (thin dashes). The 
kink becomes steeper as A is increased. 



4. Prospects 

The SSEP on a ring displays a phase transition from a homogeneous state to a kink- 
Hke profile as the trajectories' activity is lowered below its average value. We have 
described both the details of the kink profile that develops and the large deviation 
function of the activity itself in this nonuniform dynamical state. In the course of 
our analysis, we have identified three distinct finite-size scaling regimes. The system 
dealt with in the present work exhibits a second order continuous phase transition. 
It would certainly be quite desirable to find similar solvable examples for first-order 
transitions. The mechanism by which the symmetry breaking sets in may be more 
subtle. This is work in progress. 
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